/*
    <one line to give the library's name and an idea of what it does.>
    Copyright (C) 2013  Guillaume <email>

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/


#include "msSampler.h"


namespace atomism
{

bool msSampler::isSamplerRegisteredInPython=0;
msRegistrar msSampler::Registrar("msSampler", msSampler::createInstance);
 
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

std::ostream& operator<<(std::ostream& out,const msSampler& sampler)
{
 sampler.print(out);
 return out;
 }
 
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

void msSampler::registryInPython()
{
  #if USE_PYTHON
  msRegister::registryInPython();
  
  if( ! msSampler::isSamplerRegisteredInPython ) {
    
    using namespace boost::python;
    class_<msSampler::msSamplingData,boost::shared_ptr<msSampler::msSamplingData> >
         ("msSamplingData","store sampling statistics",no_init
	  );

    class_<msSampler,bases<msTreeMapper>,boost::shared_ptr<msSampler> >
         ("msSampler",
	   "(Virtual) base class for sampler",no_init
	  )
        .def("sample", (msSampler::msSamplingData (msSampler::*)( msScalarFunction&) ) &msSampler::sample,
	     "sample the function providen as argument, return a msSamplingData structure with the units of the scalar function. arg2: scalar function to sample.")
	.def("sample", (msSampler::msSamplingData (msSampler::*)( msScalarFunction& fct , const msUnitsManager&)) &msSampler::sample,
	     "sample the function providen as argument, return a msSamplingData structure using specified units. arg2: scalar function to sample, arg3: units desired.")
        .def("getNextPoint",  (bool (msSampler::*)( boost::shared_ptr<msGeneralizedCoordinates>) ) &msSampler::getNextPoint,
	     "set the next point of generalized coordinates after the current one (implementation depends on the derived class). arg2: generalized coordinates")
        .def("getRndPoint",  &msSampler::getRndPoint,
	     "set the the coordinates to a random point. arg2: generalized coordinates")
	.def("scalarProduct",   &msSampler::scalarProduct,
	     "compute the scalar product between two scalar functions. arg2: scalar function 1, arg3: scalar function 2")
	.def("begin",   &msSampler::begin,
	     "Initialize the sampler to begin the configurational scan")	
	.def(self_ns::str(self_ns::self));
      ;
    msSampler::isSamplerRegisteredInPython=1;
   }
#endif
 }

//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> msSampler::setCoordinates(boost::shared_ptr<msGeneralizedCoordinates> coors) {

  if(!coors->hasParent()) coors->setAffiliation(mySharedPtr());
  Coordinates=coors;

  return mySharedPtr();
}

//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

std::ostream& msSampler::print(std::ostream& out) const
{
 msTreeMapper::print(out);
 output::printHeader(out,"Sampler");

 if( Coordinates!=boost::shared_ptr<msGeneralizedCoordinates>() )
   {
    out<<"\nId of the generalized coordinates:  "<<Coordinates->getId()<<"; type: "<<Coordinates->getType()<<endl<<" size="<<Coordinates->noOfDim()<<"\t";
    Coordinates->getUnits()->abstract(out);
    out<<"\n\tList of unfreezed coordinates:\n";
    msGeneralizedCoordinates::UnfreezIt qi = Coordinates->beginUnfreezed();
    for(; qi!=Coordinates->endUnfreezed(); qi++)
       {
        out<<"\t\t"; (*qi).print(out); out<<endl;
        }
    }
 else out<<" no coordinates defined"<<endl;
 return out;
 }
 
 

msSampler::msSamplingData msSampler::sample(msScalarFunction& fct) {

 LOGGER_ENTER_FUNCTION_DBG("double msSamplerLinear::sample(boost::shared_ptr<msScalarFunction> fct)",
		        getFullId());

 Coordinates=fct.getCoordinates();
 AffStep = 10;
 begin();
 int width=10;
 std::stringstream comment;
 double result=0;double v=0;

 msSampler::msSamplingData data;
 Coordinates->setMinValues();

 msGeneralizedCoordinates::iterator it=Coordinates->begin();

 for(;it!=Coordinates->end();++it) comment<<std::setw(width)<<(*it)->getId();
 comment<<std::setw(30)<<"field value"<<std::setw(30)<<"integral"<<endl;
 
 LOGGER_WRITE(msLogger::INFO,comment.str());

 while( getNextPoint(Coordinates) )
      {
       result += v = fct.evaluate(); 

       if( v<data.Min ){ data.Min=v; }
       if( v>data.Max ){ data.Max=v; }
       
       if( ( nScanned!=0) && (nScanned % AffStep) == 0 ){   
          
          comment.str("");
          for(it=Coordinates->begin();it!=Coordinates->end();++it) comment<<std::setw(width)<<(*it)->Value;
	  comment<<std::setw(30)<<v<<std::setw(30)<<result<<endl;

	  LOGGER_WRITE(msLogger::INFO,comment.str());
       	  }  
       }
 data.Average=result/nScanned;
 result*=Coordinates->getVolumeDq(); 
 data.Integral=result;
 LOGGER_WRITE(msLogger::INFO,"Sampling done: "+data.abstract());
 LOGGER_EXIT_FUNCTION();
 return data;
 }
 

double msSampler::scalarProduct(msScalarFunction& fct1 ,msScalarFunction& fct2) {

 LOGGER_ENTER_FUNCTION_DBG("double msSampler::scalarProduct(msScalarFunction& fct1 ,msScalarFunction& fct2)",
		        getFullId());

 if(fct1.getCoordinates() !=fct2.getCoordinates())
    BOOST_THROW_EXCEPTION(msError("The two scalar functions must share the same coordinates",
			          "double msSampler::scalarProduct(msScalarFunction& fct1 ,msScalarFunction& fct2)",getFullId())
                          );
 Coordinates=fct1.getCoordinates();
 int width=10;
 std::stringstream comment;
 double result=0;double v=0;

 msSampler::msSamplingData data;
 Coordinates->setMinValues();

 msGeneralizedCoordinates::iterator it=Coordinates->begin();

 for(;it!=Coordinates->end();++it) comment<<std::setw(width)<<(*it)->getId();
 comment<<std::setw(30)<<"field value 1"<<std::setw(30)<<"field value 2"<<endl;
 LOGGER_WRITE(msLogger::INFO,comment.str());
 
 AffStep=10;
 begin();
 while( getNextPoint(Coordinates) ){
   
       double Fct1=fct1.evaluate();
       double Fct2=fct2.evaluate();
       
       result += Fct1*Fct2; 
       if( ( nScanned!=0) && (nScanned % AffStep) == 0 ){   
          
          comment.str("");
	  for(it=Coordinates->begin();it!=Coordinates->end();++it) comment<<std::setw(width)<<(*it)->Value;
          comment<<std::setw(30)<<Fct1<<std::setw(30)<<Fct2<<std::setw(30)<<result<<endl;

	  LOGGER_WRITE(msLogger::INFO,comment.str());
       	  } 
       }
 LOGGER_EXIT_FUNCTION();
 return result;
 }
}
